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Abstract. We consider the stability properties of the ocean of accreting magnetic neutron stars. It turns out 
that the ocean is always unstable due to the combined influence of the temperature and chemical composition 
gradients along the surface and of the Hall effect. Both the oscillatory and non-oscillatory modes can be unstable in 
accreting stars. The oscillatory instability grows on a short timescale ~ 0.1 — 10 s depending on the lengthscale of 
a surface inhomogeneity and the wavelength of perturbations. The instability of non-oscillatory modes is typically 
much slower and can develop on a timescale of hours or days. Instability generates a weak turbulence that can 
be responsible for mixing between the surface and deep ocean layers and for spreading the accreted material over 
the stellar surface. Spectral features of heavy elements can be detected in the atmospheres of accreting stars due 
to mixing, and these features should be different in neutron stars with both stable and unstable burning. Motions 
caused by instability can also be the reason for slow variations in the luminosity. 
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1. Introduction 

The surface layers of accreting neutron stars are in a 
melted state because of the heat production due to nu- 
clear burning of the accreted material. The way a neutron 
star burns the accreted hydrogen and helium to heavier 
elements is sensitive to the accretion rate (e.g., Bildstcn 
1998). This nuclear burning is very to be stable for rel- 
atively high accretion rate, M > 2.6 x 10 _8 M Q yr _1 . 
Burning at such accretion rates is not hot enough to pro- 
duce iron group or heavier elements; thus the surface lay- 
ers of neutron stars consist mainly of CNO group elements 
(Brown & Bildsten 1998). Fusion of this fuel to heavier el- 
ements does not occur at the density p < 10 9 g cm -3 . If 
the neutron star is accreting pure He then the burning at 
p > 10 9 g cm~ 3 produces iron group elements, but if the 
star accretes a mixture of H and He then the products are 
much havier (Schatz et al 1999). The low ionic charge and 
a relatively high temperature (> 10 8 K) caused by the nu- 
clear burning delay cristallization until very high densities 
are reached. Therefore neutron stars accreting at a high 
accretion rate are covered with massive oceans. Neutron 
stars accreting at lower rates, M < 2.6 x lO -8 M0yr -1 , 
burn accreted H and He directly onto iron group elements 
(Lewin, van Paradijs & Taam 1992, Schatz et al. 2001, 
Woosley et al. 2004), while the Coulomb interaction be- 
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tween heavy ions leads to crystallization at much lower 
densities and to a much shallow ocean. 

MHD processes in the ocean are of general interest 
to studies of magnetic and thermal evolution, thermonu- 
clear burning, neutron star seismology, and properties of 
quasi-periodic oscillations. For instance, in weakly magne- 
tized accreting neutron stars, the ocean of light elements 
supports shallow water waves, which can be the reason 
for low frequency (~ 5 — 7 Hz) quasi-periodic oscillations 
in the brightest X-ray sources (Bildsten & Cutler 1995, 
Bildsten, Ushomirski & Cutler 1996). Under certain con- 
ditions, large-scale flows in the ocean can influence the 
magnetic evolution of accreting neutron stars (Bisnovatyi- 
Kogan & Romberg 1974, Rai Choudhuri & Konar 2002). 
MHD phenomena can play a particular role during ther- 
monuclear X-ray bursts. Thus, if an unstable nuclear burn- 
ing is not spherically symmetric, then burning fronts can 
spread around the star by igniting the accumulated cold 
fuel ahead of the moving front (Fryxell & Woosley 1982, 
Bildsten 1993, Spitkovsky, Levin & Ushomirsky 2001). 
The propagation speed can be considerably increased 
when convection and enhanced mixing are occurring at 
the burning front since convective motions are very ef- 
ficient at transporting heat in the ocean. Note that ef- 
fects caused by rotation should be important during X- 
ray bursts since most bursters are rapidly rotating. Most 
likely, an expanding ocean rotates differentially, and dif- 
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ferential rotation can be the reason for MHD instabilities 
(Cumming & Bildsten 2000, Menou 2004). 

Mixing in stars is often atributed to hydrodynamic in- 
stabilities, such as convection. Standard thermal convec- 
tion arises if the temperature gradient is superadiabatic 
and occurs only in the atmosphere of hot NSs with a 
weak magnetic field (Miralles, Urpin & Van Riper 1997). 
Convection, however, is not the only instability that oc- 
curs in the ocean. Stability properties of accreting neu- 
tron stars can be complicated because of the presence of 
a strong magnetic field, B ~ 10 12 — 10 13 G, and of the 
thermal and compositional gradients that generally are 
not parallel to gravity g. A possibile occurrence of the 
Kruskal-Schwarzschild instability in polar caps of accret- 
ing neutron stars has been considered by Litwin, Brown & 
Rosner (2001), who found that ballooning modes with the 
displacement perpendicular to the magnetic field can be 
unstable if the overpressure at the bottom of the neutron 
star ocean exceeds magnetic pressure by a factor ~ 8a/ h, 
where a and h are the horizontal lengthscalc of the cap and 
the density scale hight, respectively. The crucial point of 
the instability is line-tying to the neutron star crust and, 
therefore, unstable modes are localized to within a density 
scale height of the ocean bottom. The ballooning instabil- 
ity can arise on a timescale ~ 30 hr if M <~ 10 _10 M Q /yr 
and, most likely does not affect the nuclear burning but 
can degrade the confinement to prevent accumulation of 
mass above the instability threshold. 

One more instability that can manifest itself in ac- 
creting neutron stars is the well-known Parker instabil- 
ity (Parker 1966). The downward flow of accreted matter 
may generate strong non-uniformity in the magnetic field 
and, hence, the electric current, which in turn can lead 
to Parker-type instability (Cumming, Zweibel & Bildsten 
2001). Note that this instability arises only if both the 
magnetic field and its non-uniformities arc sufficiently 
strong and if the field decreases with height. Employing 
the effect of thermal diffusivity, Cumming, Zweibel & 
Bildsten (2001) estimated that the critical field above 
which the instability occurs is ~ 10 12 G at the top of 
the crust. The Parker- type instability can be of particular 
importance if the magnetic field in the surface layers has 
a complex geometry, for example, such as one considered 
by Urpin & Gil (2004). 

The material accreted from a companion star car- 
ries angular momentum. The inflow of angular momen- 
tum promotes differential rotation in the ocean of neu- 
tron stars and can cause rotation induced instabilities. 
These instabilities occur when the angular velocity gra- 
dient is so steep that the destabilizing effect of shear over- 
whelms the stabilizing effect of buoyancy (see, e.g., Zahn 
1983). In the context of neutron star oceans, the instabil- 
ities caused by differential rotation have been considered 
by Fujimoto (1988, 1993). Turbulent motions generated 
by the instabilities transport the angular momentum in 
the ocean effectively, but the elemental mixing is much 
less efficient because it is hampered by a stable stratifica- 
tion. Note that a stable stratification can suppress some 



rotation induced instabilities since the Brunt- Vaisala fre- 
quency is typically much larger than the angular velocity 
in accreting neutron stars (Cumming & Bildsten 2000). If 
differential rotation departs from the cylindrical symme- 
try, then the ocean can be subject to baroclinic instability 
(Cumming & Bildsten 2000). A stable stratification gen- 
erally suppresses the modes with horizontal lengthscales 
that are much longer than a pressure scale height, but 
short-wavelength can be unstable, and turbulent trans- 
port caused by the baroclinic instability can be significant 
during type 1 X-ray bursts. Note also that some mixing 
in accreting neutron stars can be provided by tidal waves 
(Lou 2001). It is unlikely, however, that this mechanism is 
efficient since tidal waves in the neutron star ocean shoud 
have a small amplitude because of high gravity. 

In the present paper, we consider the instability caused 
by horizontal inhomogeneities of the temperature and 
chemical composition in the ocean of accreting magnetic 
neutron stars. The thermal and compositional state of 
the ocean departs substantially from the spherical sym- 
metry for such stars. The departure can be caused, for 
example, by the accreted hydrogen and helium focusing 
onto the magnetic poles. It is usually assumed that the 
accretion flow interacts with the magnetosphere at the so- 
called Alfven radius, Ra, which is determined by the bal- 
ance of the magnetic and inertial forces. For the standard 
neutron star with a dipole magnetic field and with mass 
M = 1.4M Q and radius R w 10 6 cm, the Alfven radius is 
given by 

i?A~2x itfB^MZW cm : (!) 

where B 13 = B/10 13 G and M_i = M/(10- 10 M Q yr- 1 ). 
For a weakly magnetized star accreting at a high accre- 
tion rate, the Alfven radius becomes comparable to the 
stellar radius, and the magnetic field cannot funnel the 
accreted material onto the magnetic poles. Therefore, the 
distribution of the accreted material over the stellar sur- 
face is more or less smoothed in such stars. In contrast, 
a strongly magnetized neutron star accretes the mate- 
rial predominantly onto the magnetic poles, and distribu- 
tion of the chemical composition and temperature is sub- 
stantially non-spherical in the ocean of magnetized stars 
(Brown & Bildsten 1998). 

It has been shown by Urpin (2004) that the combined 
influence of the Hall effect and a thermal gradient perpen- 
dicular to g can be responsible for instability in the ocean 
of isolated neutron stars. In accreting neutron stars, the 
presence of a compositional gradient that generally is also 
non-parallel to g can lead to similar instability. This sort 
of instability is qualitatively different from instabilities al- 
ready considered in accreting neutron stars because it oc- 
curs only in the presence of the Hall effect. The criterion of 
instability can be satisfied in many accreting neutron stars 
if the compositional gradient is not parallel to gravity. 

This paper is organized as follows. In Section 2, we 
consider the main equations governing hydrodynamic mo- 
tions in the ocean of accreting magnetic neutron stars and 
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derive the dispersion equation. Section 3 deals with the 
criterion and growth time of instability. In Section 4, we 
discuss mixing processes associated to compositional in- 
stability in the ocean. 

2. The basic equations 

We consider the linear stability of the ocean on an accret- 
ing neutron star assuming that both the temperature and 
chemical composition can generally depart from spherical 
symmetry in the unperturbed state. We do not specify the 
mechanisms responsible for these departures but they can 
originate, for example, in the magnetic field B that funnels 
the accretion flow onto the poles and provides anisotropy 
to the heat transport. As a result, gradients of the temper- 
ature and chemical composition are not parallel to gravity 
g, which is approximately radial in the ocean. 

The depth of the ocean depends crucially on the tem- 
perature and chemical composition and varies over a wide 
range. Crystallization of a Coulomb liquid occurs when 
the ion coupling parameter T = Z 2 e 2 / (aksT) reaches the 
critical value r = T m w 170 - 180 (Slattery, Doolen & De 
Witt 1980, Potekhin & Chabrier 2000); a = (3/4 7 m J ) 1 /3 
is the mean inter-ion distance, n, and Z are the number 
density and the charge number of ions, respectively; ks 
is the Boltzmann constant. Then, the crystallization tem- 
perature is 



Z 2 e 2 
ak B T m 



1.3 x 10 5 Z 5 /V/ 3 (—) K 



(2) 



where x = Zp/10 6 A g/cm 3 , and A is the atomic number 
of ions. In neutron stars with temperature T > 10 8 which 
is more or less typical for stars accreting with the rate 
<~ 10 _9 M Q yr -1 , crystallization occurs at the density p <~ 
3 x 10 10 g/cm 3 if Z = 8. 

The equations governing velocity, magnetic field, and 
thermal balance are (see, e.g., Landau & Lifshitz 1959) 



v + (v ■ V)v 
p + V • (pv) = 0, 



Vp 1 

— +g + —{VxB)xB, 

P iTTp 



B - V x (v x B) = -V x [i) ■ (V x B)}, 
V • B = 0, 

t - —p + v • (AVT) = —V • (k ■ VT), 
pc p pc p 



(3) 

(4) 
(5) 
(6) 

(7) 



where v is the fluid velocity; AVT = VT - V a dT is the 
difference between the real and adiabatic temperature gra- 



thermal expansion coefficient and c p the specific heat at 
constant pressure, p; and finally fj = c 2 i?/47r, R and k are 
tensors of the electrical resistivity and thermal conductiv- 
ity. Following some tensor operations on Eqs. (5) and (7), 
we have 

k- VT = k (0) 6(6 • VT) + k (b) [VT - 6(6 • VT)] + k (a) 6 x VT, 



where and are the tensor components along and 
across the magnetic field, respectively, is the so called 
Hall component, and 6 = B/B. An analogous expression 
can be written for fj. Note that for the magnetic diffu- 



sivity we have rj^ 
parameter, 



r^( B ) = ?/ a )/q! where a is the Hall 



a 



n 



BeT 



9.9 x 10 3 Bi 3 
ZAil + x 2 / 3 )' 



(8) 



Qbc and t are the electron gyrofrequency and relaxation 
time, respectively; and A is the Coulomb logarithm. 

We neglect the viscous term in the momentum equa- 
tion (3) since viscosity plays a less important role among 
kinetic processes in the ocean. Usually, both electrons and 
ions contribute to the shear viscosity. The ion shear viscos- 
ity dominates at a relatively low density but is typically 
much smaller than magnetic diffusivity r\ and thermal dif- 
fusivity, x — k /p c p- Electron viscosity is greater in deep 
layers of the ocean (Itoh, Kohyama & Takeuchi 1987). 
The ratio of electron viscosity v e and magnetic diffusivity 
along the magnetic field is given by 



n 



v e _ 2.2 x 10 2 



,5/3 



AZA (1 + x 2 / 3 ) 2 ' 



(9) 



where A is the Coulomb logarithm. The kinematic vis- 
cosity becomes comparable to the magnetic one only at 
a high density > 10 9 g/cm 3 if the ocean consists of light 
elements. Therefore, neglecting viscosity in Eq. (3) seems 
to be qualitatively justified. 

We assume that the ocean is in hydrostatic equilibrium 
in the unperturbed state, 



— = G = g + ——(V x B) x B. 

p 47T/9 

Taking the curl of this equation, we have 
X7px G = --5-V x [(V x B) x B}. 



(10) 



(11) 



The density and, hence, temperature and compositional 
gradients have components perpendicular to G if the mag- 
netic field is not force-free in the ocean. 

Consider only relatively fast hydrodynamic processes 
varying on a timescale shorter than the characteristic 
timescale of nuclear reactions. The time to burn all hydro- 
gen into CNO-group elements is, for example, on the order 
of 20 hours at high accretion rates (see, e.g., Cumming & 
Bildstcn 2000). Under this assumption, the number den- 
sity of a species j satisfies the continuity equation 



dicnts, V ad T = /3Vp/pc p ; j3 = -{d\n p/ d\nT) p is the Q n 



-^+V-(n j «)=0. 



(12) 



For our purposes, it will be convenient to characterize the 
composition by a mass fraction of species Xj = pj/p, 
where pj is the density of species j. Note that if the num- 
ber of species is N then only N—l fractions Xj are linearly 
independent, since Xj satisfy the condition ^2 N Xj = 1; 
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summation is over all species of ions. If the chemical com- 
position is frozen in a moving fluid parcel, then the equa- 
tions governing Xj read 



9X± 
dt 



+ V • (Xjv) = 0. 



(13) 



The equations for small perturbations can be ob- 
tained by linearization of Eqs. (3)-(7), and (13). We as- 
sume that in the unperturbed state, a steady tempera- 
ture and composition gradients are maintained, and there 
are no motions. In what follows, small perturbations will 
be marked by the subscript 1, but subscripts will be 
omitted for unperturbed quantities. We consider short- 
wavelength perturbations with spatial and temporal de- 
pendence exp(7i — ik ■ r) where k is the wave-vector. 
Linearization of Eqs. (3)-(7), (13) yields 



1Vl = z fc^l + G— - — (k x Br) x B, 
p p inp 

7Pi + k ■ v\ = 0, 



(14) 
(15) 



(7 + w,)Bi = -i(k ■ B)v 1 - 77a (fc • b)k x B u (16) 
fe-Bi=0, (17) 



7/3 

(7 + oj k )T x - -iz-pi = -vx • (AVT), 
pc p 

7X1J + Vl ■ VXj = 0, 



(18) 



(19) 



where the characteristic inverse timescales of the ohmic 
dissipation and thermal diffusion are given by u; n — r/^k 2 
and u K = [n^k 2 + (k^ - «( B ))(fe • b) 2 ]/pc p , respectively. 

Perturbations of the pressure in Eq. (14) can be ex- 
pressed in terms of perturbations of the density, temper- 
ature, and fractions of species X\j, using the equation of 
state, 




N-l n 



where Sj = — (din p/dXj) P T are the coefficients of chem- 
ical expansion, and summation is over N — 1 linearly in- 
dependent species. 

Consider Eqs. (14)-(20) in the case kc s 3> 7 that corre- 
sponds to the Boussinesq approximation for slowly varying 
modes. Then, the dispersion equation reads 



7 + a 5 7 5 + 047 + a 3 Y + a 2 7 + ai7 + a = 



(21) 



where 
a 5 = 

(24 = 



uj 2 +u 2 a + 2u rj u K + 2oj a — u , 



a 3 = lu k (uj 2 +lu 2 a + 2uj 2 a - uj x ) + 2uj ri (uJ 2 A - lo 2 ), 

o I 2 2 \ 2^2 2 1 2\, 4 

fl 2 = 2lV n LJ K (uJ A - UJ X ) - W K, + UJ A + UJ A ) + UJ A , 

2 [ 2 1 2 , 2 1 2 / 2 1 2 1 2 \ 

a x = uj a [uj k uj a - ijJ v [u! + uj h \ - uJ K uJ Y (uj, q + uj a + oj a ), 
a = -uj k uj v uj a (uj x + u) 2 XH ). 



The characteristic frequencies in this equation are 

(k ■ B) _ 
w A = arj\\k(k ■ b) , oja — — 7^=" i = D ' 



uj 2 = ■ AVT , uj 2 =uj 2 g +u 2 x =pD-C, 

"H = «%H + O&H > = ^ AVT ■ (fe X D) 



k 2 T 



W L = ^(fexD).VX, 



where vectors D and C are given by 



D = G- ¥ {k-G), C = — + — , 

and we denote VX = VJ^ -1 SjVXj. In the limit ~ 
wxh w 0, Eq. (21) yields the dispersion equation for a 
chemically homogeneous ocean derived by Urpin (2004). 

3. The growth rate and criterion of instability 

Equation (21) describes six modes that exist in a chem- 
ically inhomogeneous ocean of magnetic neutron stars. 
Stability properties of the modes depend crucially on k 
and, under certain conditions, modes can be unstable. The 
dissipative frequences w v , uja, and uj k are typically smaller 
than the dynamical frequences oja, ojo, and ojh except 
perturbations with a very short wavelength, A = 2-ir/k. 
Therefore, the dispersion relations for some modes can be 
obtained by making use of the perturbation method and 
assuming that dissipative effects are small. We expand the 
growth rate as 7 = 7^ + + + ... where 7^, 7^, 
and 7 1 - 2 -' are terms of the zeroth, first, and second order 
in dissipative frequencies, respectively. The corresponding 
expansion should be made for the coefficients of Eq. (21) 
as well. In the zeroth order when dissipation is neglected, 
Eq. (21) reduces to a quadratic equation, 



(20) 7 (0) 4 + (2-i - ^)7 (0) 2 + - <4) = 0, 



(22) 



This equation describes four modes with the frequencies 

(23) 



(0) 2 2 (0) 2 2,2 

7l,2 =-WA, 73.4 =-U A +U Q . 



The fifth and sixth roots of Eq. (21) are vanishing at 
the zeroth order of approximation, 7^ = = 0. 

Usually, 7^°) 2 < for the modes 1-4 because the inequal- 
ity uj 2 > uo 2 A requires temperature and composition gradi- 
ents that do not exist in neutron stars. Therefore, modes 
1-4 are oscillatory, and their instability does not occur at 
the zeroth order. 

For applicability of the perturbation procedure, the 
first corrections to 7^ caused by dissipative effects have 
to satisfy the condition |7^| <C |7^| ~ wa, since uja 
is usually larger than other frequencies. The real parts of 
for the modes 1-4 have a particularly simple form if 
I I and \u)h\ satisfy the inequalities 



\lo 2 \ > lo v lo k (1 + a 2 ) , \u>%\ > lo v uj k (1 + a 2 ). 



(24) 
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Under this condition, the real parts of 7W are given by 



Re 7$ 
Re 7 $ 



ri (LO 2 H 



2ul 



u v u a [uj h + uj (1 — C)] 



2w 2 (o 



(25) 
(26) 



where £ = w K Wj/w,w^. Note that if condition (24) is not 
fulfilled and dissipative effects dominate the dynamical 
influence of temperature and composition gradients along 
the surface, then modes 1-4 are stable. 

The roots 5 and 6 of Eq. (21) describe a couple of non- 
oscillatory modes and arc linear in dissipative effects. To 
calculate these roots we should keep the terms in Eq. (21) 
up to the second order in dissipative frequencies assuming 
that 7 is linear at these frequencies. Then, roots 5 and 6 
are approximately 



75,6 ~ 

where 



A 



± 



4 



A , 



(27) 



LO K LO. q (uJ 



X 



UJ K {UJ 2 A - LU 2 X ) - U V (U)1 + 



A 1 = 

Modes 5 and 6 are secular and appear only due to dissipa- 
tive effects. Note that contrary to modes 1-4, expression 
(27) for non-oscillatory modes applies even for small wave- 
lengths which do not satisfy the inequality (24). The only 
condition for applying Eq. (27) is <C uja- 



3.1. Instability of the oscillatory modes 

The real parts of oscillatory modes are given by Eqs. (25) 
and (26). Note that the vertical component of VX can 
be rather large because the chemical composition varies 
substantially with the ocean depth. Therefore, the contri- 
bution of VX to ujq can usually be comparable to (or even 
larger than) that of AVT. As mentioned, usually lo\ > u>q 
(see Eq. (23)). If Uq > then instability of oscillatory 
modes occurs if 



uj 2 h > uj 2 or uj% < uj 2 (C - 1). 



(28) 



In contrast, if stratification is convectively stable and uj 2 < 
then these modes are unstable if 



uj 2 h < uj 2 or uj% > uj 2 (C - 1). 



(29) 



In the particular case when gradients of the temper- 
ature and chemical composition arc parallel to G and 
u) 2 H = lo 2 xh = 0, modes 1 and 2 are stable but modes 
3 and 4 can be unstable if £ > 1, or 

oj k lj 2 > lu v u 2 a . (30) 

This condition generalizes the well-known criterion of os- 
cillatory convection (see, e.g., Chandrasekhar 1961) for the 



case of a chemically inhomogencous fluid. In the ocean of 
neutron stars, parameter £ can be estimated as 



C ~ 3 x 10" 



T 7 \ 2 x 2 [l + {l + x 2 / 3 ) 1 ' 2 ] 



(31) 



ZA 2 B 2 3 H 3 (l + .T 2 /3) 2 

where H = C^ 1 = {H T X + H^ 1 )- 1 , H T = \(AVT/T)^\, 

and Hx = /3|(VX)ii l^ 1 are the lengthscales of the temper- 
ature and chemical composition along G (the subscript | 
denotes the component parallel to G), H3 = H/10 3 cm; 
and A2 = A/ 100 cm. Usually £ < 1, and condition (30) is 
not satisfied except for very hot stars with a low magnetic 
field (Miralles, Urpin & Van Riper 1997). 

If (CI < 1, then conditions (28) and (29) read 



oj 2 h |>| tul 



(32) 



This inequality can always be satisfied by a corresponding 
choice of k once (AVT)i ^ or (VI)i ^ 0, where _L 
denotes the component perpendicular to G. Introducing 
local coordinates with the z-axis antiparallel to G and the 
a;-axis aligned with k±, we have G — —Ge z , k = k x e x + 
k z e z ; e x ,e y , and e z are unit vectors. The expressions for 
w 2 and ujjj read in these coordinates 

lu 2 = j^f3G(k z C x - k x C z ) , uj 2 h = ^a(3G(k ■ b)C y , (33) 

where C X:VjZ are the corresponding components of C. 
Then, condition (32) yields 



a I k ■ b 



Cy I > I k x C z k z C x 



(34) 



The instability occurs only if k has a non-vanishing com- 
ponent perpendicular to the plane (G,C). Perturbations 
are suppressed if k is parallel to this plane and, hence, 



Cy=0. 



Equation (34) applies only if condition (24) is satisfied 



and 

I k x C z — k z C x I 3> 



l3Gk x 



(1 + a 2 ). 



(35) 



If inequality (35) is satisfied then we have for unstable 
perturbations 

a\k-b\C v »^^(l + a 2 ). (36) 

Note the important difference between Eqs. (34) and (35). 
Condition (34) is linear in k, so it can always be satisfied 
by proper choice of the direction of k. On the contrary, 
inequality (35) is non-linear and yields restrictions on the 
wavelength of unstable perturbations rather than on their 
direction. 

Consider initially condition (34). This condition is ob- 
viously fulfilled if the wavevector is close to the direction 
determined by 

k x C x H 
k z C z L 

where L = CJ 1 - (i T x + i^ 1 )" 1 , L T = |(AVT/T) J^l, 
and Lx = /3|(VF)j_| _1 are the temperature and compo- 
sition lengthscales perpendicular to G. One can estimate 



(37) 
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the angle AO around this direction where inequality (34) 
is still satisfied. Restricting ourselves by linear terms in 
AO, we have from Eq. (34) 



b z + b x 



Note that the components C x and C y are usually small 
compared to C z , 



I a 



^ » 10,1 



\CJ ~ - 



Therefore, Eq. (38) is approximately equivalent to 
ICJ ,. ,H 



AO < a\b z 



\C Z 



a\b z 



(39) 

(40) 
w 0. 



except the region near the magnetic equator where b z 
This estimate is correct if a < L/H, or 

B < 2 x 10 n ZL 5 Hs 1 (l + x 2 ' 3 ) G, (41) 

where L5 = L/10 5 cm. Note, that for all perturbations 
with k within angle AO we have k x < k z . If magnetization 
is very strong and a > L/H then inequality (34) can be 
satisfied even for perturbations with k x ~ k z . 

Condition (35) determines the wavelength of unstable 
perturbations. Estimating C x <~ 1/L, we can transform 
Eq. (35) into 

1 



> 



f(a), f(a) = l + a 2 



(42) 



L ' k x f3G 

(we assume that k is not perpendicular to B). Using the 
definitions of oj v and uj k and taking account of Eq. (37), 
this expression can be rewritten as the condition for the 
wavelength, 



*>A--^(£n^)"> ( <. 



0.6/ 



1/4 



Z 2 L 2 



1/4 



[l + a + x 2 / 3 ) 1 ^ 2 ]- 1 / 4 



cm; 



(43) 



\H 3 x 1 / 3 / 

we use the electron thermal conductivity by Urpin & 
Yakovlev (1980). For applicability of a short wavelength 
approximation, A and A cr should be smaller than the verti- 
cal lengthscale in the ocean, H <~ 10 — 20m. This condition 
can in general be satisfied for a wide range of parameters 
except in the case of a very strong magnetic field. In such 
field a ^> 1, and we have / w a. Then, critical wavelength 
X cr becomes larger than H if 



B > B d ~ 4 x 10 1 



tt \ 1/2 

fr) x 1/6 d + . 2/3 ) 
x[1 + (1 + .t 2 / 3 ) 1 / 2 ] 1 / 2 G. 



(44) 



Critical field Bd is much stronger than typical magnetic 
fields of the majority of pulsars. 

The growth rate of instability depends on the wave- 
length and can vary within a wide range. If lo\ » uj 2 and 
C < 1, then Re7 for unstable modes is given by 

a(kb)C y 



k x C z k z C x 



.(45) 



The growth rate is maximal if k lies in the plane perpen- 
dicular to the plane (G, C). Then C x = and, estimating 
C y /C z ~ H/L, from Eq. (45) we have 



(38) R c'- mk2a ^ Cy 



k x C z 



IQ-^B^e^x 



(46) 



The growth rate increases rapidly with increasing k and 
reaches its maximum at A <~ A cr . According to Eq. (37), 
the minimal value of the ratio e = k x /k is on the order of 
H/L. Substituting this estimate into Eq. (46) and assum- 
ing, for example, that A2 ~ 1, we have 



Re7 



0.1 Bigx" 1 s _1 . 



(47) 



The growth time for such perturbations is ~ 10 s in the 
layer with p ~ 10 6 g/cm 3 if B ~ 10 13 G. Perturbations 
with a shorter vertical wavelength can grow even faster. 
For example, perturbations with A ~ A cr ~ 10 cm grow 
on a timescale of about 0.1 s. 

3.2. Instability of the non-oscillatory modes 

Root (27) describes a couple of non-oscillatory modes that 
can generally be unstable in the ocean. The instability 
conditions of these modes are 



A < or A 1 < 0. 



(48) 



If we assume oj\^> uj 2 , which is typical for those neutron 
stars with B > 3 x 10 s x 1 / 2 X 2 H^ /2 L^ 1 G, then conditions 
(48) read 



u 2 x + 4fl > 0, 



or 



u k (lu 2 a 



< 0. 



(49) 



(50) 



Consider initially the condition (49). If the magnetic 
field is weak or G and VX are almost parallel (or antipar- 
allel), we have \u)x\ ^> \uj"x H \, and inequality (49) trans- 
forms into lo 2 x > 0. This condition has exactly the same 
form as the standard criterion of the Ledoux convection. 
Note, however, the big difference because the standard 
Ledoux convection is a dynamical instability, whereas the 
instability considered in this subsection is secular. Using 
the local cartesian coordinates, we can represent condition 
uj 2 v > as 



N ~ 1 rlX 

V..V »'^< 



dz 



(51) 



As an example, in the simplest case of a mixture of two 
species with the atomic masses mi and m 2 , Eq. (52) reads 
8\dX\/dz < 0. If m2 > tni then 5\ = —{dp/dXi) P T is 
positive, and the condition of instability is dXi/dz < 0, or 
d{n\/n2)/dz < where n\ and ni are the number density 
of species 1 and 2, respectively. Therefore, the instability 
occurs if the relative number density of a light species 
increases with depth. Most likely, however, the relative 
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concentration of heavy elements increases with depth and 
uj\ < 0. Therefore, the instability most likely does not 
occur in the ocean if VX || G. 

If VX and G are not parallel, then the stability prop- 
erties change drastically. In the general case, inequality 
(49) yields 

G VX - p (fc ■ G) (fc ■ VX) - a( ^ 2 b) G-(kxVX) > 0.(52) 

Again using the local cartesian coordinates, we have from 
this equation 



k x k z (ab z V y X + V X X) > k 2 x {V Z X - ab x V y X), 



(53) 



where V x ,y,zX are cartesian components of VX. If the 
composition varies over the surface, then for any distri- 
bution of chemicals, there exists a region in the (k Xl k z )- 
plane where the condition of instability (53) is satisfied. 
The non-oscillatory modes become unstable if even one 
of the components V X X or V y X is non-vanishing. If the 
vertical component of a chemical gradient is large and 
V Z X > ab x V y X, then the instability occurs if 



< 



ab z V y X 



V Z X - ab x V y X 



(54) 



If V Z X < ab x V y X, then the condition of instability is 



> 



ab z V y X 



v x x 



v z x 



ab x V y X 



(55) 



Note that the instability occurs even if the magnetic field 
is weak and the magnetization parameter a is small, a <C 
1. In this case, a surface chemical inhomogeneity can be 
caused, for example, by accretion from the disc and may 
have a belt-like structure. If a <C 1, then from Eq. (54) 
we have 



< 



V X X 



v z x 



E2L 



(56) 



and the instability arises only for perturbations with the 
vertical wavelength much shorter than the horizontal one. 
In a moderate magnetic field with Lx/Hx > a > 1, in- 
equality (55) yields 



< 



ab z VyX 



v z x 



By 



el- 



and the region of unstable wavevectors is more extended 
than where a <C 1. If a > Lx/Hx, then the condition of 
instability is given by Eq. (56) and reads 



> 



1. 



(58) 



Note, however, that the condition a > Lx/Hx requires a 
rather strong magnetic field, 



B > 2 x 10 n Z(l + x 2/3 ) 



2/3 \ L X5 



H 



G. 



A'.'i 



(59) 



which is impossible in many accreting pulsars {Hxz = 
Hx/10 3 cm, L X5 = L x /10 5 'cm). 



Condition (50) generalizes the Chandrasekhar crite- 
ria of stationary convection for the case of a chemically 
inhomogeneous fluid with non-parallel gravity and com- 
position gradient. Since u>a is typically much larger than 
other characteristic frequencies and uj k S> uj v , it is unlikely 
that inequality (50) is fulfilled in "standard" neutron stars 
because these conditions require very low magnetic field 
and large chemical gradient. 

The growth rate of the composition-driven instability 
(A < 0) is always proportional to the gradient of X along 
the surface. Taking into account that u K 3> u) v and assum- 
ing that the Alfven frequency is much greater than any 
other characteristic frequency, we obtain \Ai\ 2 ^> \Ao\. 
Then, roots 5 and 6 are approximately given by 



(i) 



A 



(i) 



-A 1 . 



(60) 



As it was mentioned, the condition A\ < is most likely 
not fulfilled in the ocean, which means mode 6 is usually 
stable. Mode 5, however, can be unstable if A a < 0. Since 
uj\ » \uj 2 \ the growth rate of this mode is 



75 



(i) 



-4-{u x +Jxh) = 



4:TrpGrj k x 
(fc • B) 2 



[a{k-B)V v X-k x V z X + k z V x X\. (61) 



We estimate the growth rate for the case when fc lies in 
the plane perpendicular to the plane (G, VX) and, hence, 
(VX) X = 0. Assuming that the magnetic field is not very 
strong {Lx/Hx > a) and the wavevector satisfies the 
condition of instability (54), we obtain from Eq. (61) 



(i) 



k x G_ 
Ifc • 61 iv. 



\VyX\ 



10" 



;A k x g 14 



pi 



Z k z BizLxh 



\ (62) 



where oj pi = ZeB / Am p c is the ion gyrofrequency and 
5i4 — ff/10 14 cm/s 2 . In contrast to the case of oscillatory 
modes, the growth rate of a non-oscillatory mode depends 
on the direction of a wavevector rather than on its value. 
Instability is faster in a neutron star where the surface in- 
homogeneities of a chemical composition are stronger. If 
the magnetic field is not very strong and satisfies the con- 
dition Lx/Hx > a > 1 (or 10 n ZA(l + x 2 / 3 )L X5 /H X5 > 



(57) B > 10 9 ZA(1 



r 2/3 



)), then the non-oscillatory mode 



grows on a relatively short timescale <~ 1 d if L ~ 100 
m. Note that the chemical and thermal gradients are very 
likely to be larger in accreting neutron stars than in iso- 
lated ones and one could thus expect more efficient hydro- 
dynamic motions in the ocean of accreting stars. 

4. Discussion 

It turns out that the ocean of neutron stars is subject to 
different instabilities if the temperature or chemical com- 
position varies over the surface. The main driving forces 
of instability are the Hall effect and horizontal advection 
of heat or composition; as a result, the growth rate of 
instability is proportional to the Hall parameter a and 
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the component of VT or VJ along the surface. This 
point can be clarified by a simple qualitative considera- 
tion. Consider the case of oscillatory modes in the magne- 
tized ocean (a > 1) with strong chemical inhomogeneities 
(|VX| > 0\&VT\/T). Then, the amplitude of unstable 
magnetic perturbations is changed mainly by the Hall ef- 
fect. For example, amplitude B\ y increases after At by 



AB 



At 



r) A (kb)(k x Bi)j, 



k 2 

VA(kb)—B lz 



(63) 



(using the divergence- free condition (6) for £?i), and small 
perturbations are marked by the subscript "1". Since os- 
cillatory perturbations are approximately Alfvenic in the 
ocean, the components of velocity and magnetic field are 
related by B Xz ~ B\ y (v\ z / V\y) . For the oscillatory in- 
stability, perturbations of the composition are relatively 
small, X\ cx (v\ ■ VX) w 0, then V\ z jv\ y ~ VX y /VY z . 
Substituting these expressions into Eq. (63), we obtain the 
growth rate (46). It is seen from this simple consideration 
that the Hall effect (Eq. (63)) is one of the driving forces of 
instability which makes qualitatively different from other 
instabilities in accreting stars. 

It is worth noting that we neglect rotation in our sta- 
bility analysis. Rapid rotation will often change the stabil- 
ity properties of accreting neutron stars. The dispersion 
equation (21) is valid until angular velocity £1 is small 
compared to all dynamical frequencies. Since ua is larger 
than loq in many cases of interest, our analysis applies 
if m.m(u>o, u>h) > ^- This condition can be written as 
<C y/a/(l + a)oj a or, using expression (33) for uj , as 



P>2xlQ- 



i L 5 



1 + a 



a 



(64) 



where P is the rotation period. 

The growth time of both the oscillatory and non- 
oscillatory instabilities is relatively short, and instabili- 
ties are likely to operate operate in a non-linear regime. 
We can estimate the saturation velocity using the mixing- 
length model (e.g., Schwarzschild 1958) that assumes that 
the turnover time of turbulence generated by instabil- 
ity is on the order of the growth time of this instability. 
Consider, for example, the saturation regime of oscillatory 
modes. Then, in a turbulent cell with the characteristic 
vertical lengthscale A z , the vertical velocity component in 
saturation can be estimated as vt z (A z ) ~ A z Re7(A z ). 
Using Eq. (46), we have 



vtz ~ 0.1B 13 x 1 e 1 A^ 2 1 — - cm/s. 

L5 



(65) 



Since the instability is anisotropic, the saturated turbu- 
lent velocity should be anisotropic as well. From Eq. (4), 
we can estimate the turbulent velocity along the surface as 
vtx ~ vt z {X x /\ z ) where \ x is the characteristic length- 
scale of turbulence along the surface. 

Turbulent motions can enhance transport in the ocean 
both vertically and horizontally. For example, turbulence 



can efficiently mix the material in the ocean. The coeffi- 
cient of turbulent diffusion can be estimated as a prod- 
uct of the turbulent lengthscale and velocity in the corre- 
sponding direction. Then, we have 



vtz ~ vtzX z ~ lOB^e x — — cm /s (66) 

L5 

for the coefficient of vertical diffusion, and 



v Tx ~ v Tx \ x ~ v Tz 



A? 



I OL'j \r ^| cm 2 /s (67) 



U \ 2 



for the coefficient of diffusion along the surface. Assuming 
that e w k x /k z <~ H/L and Bi 3 = 0.1, x = 1, we obtain 
vtz ~ 10 2 cm 2 /s. Such diffusion is sufficient for mixing the 
outer layer with depth ~ 10 m on a timescale ~hours. The 
rate of vertical turbulent diffusion should be compared to 
the gravitational sedimentation in the ocean. The coeffi- 
cient of interspecies diffusion (Brown, Bildsten & Chang 
2002) reads in our notations as 



V 



4 x 10" 



A°- 5 Z°- 7 Z$- 3 x- 



rpl.2 

1 7 2 - 

cm s 



0.6 



(68) 



where Z 2 is the charge of a trace component. This quantity 
is usually much smaller than vtz- As a result, heavy ele- 
ments can be transported by turbulence from deep ocean 
layers to the surface and manifest themselves in spectra 
of both isolated and accreting neutron stars. Spectral fea- 
tures should be different in spectra of accreting neutron 
stars with stable and unstable burning. If the accretion 
rate is high and the burning stable, then spectral features 
corresponding to CNO group elements could be detected 
in such stars. For example, these features could be ob- 
served in spectra of strongly magnetized accreting neutron 
stars where the burning is stable and a deep ocean of CNO 
group elements is formed (see, e.g., Bildsten 1998). On 
the other hand, if burning is unstable and the star burns 
H and He directly to iron group elements, then mixing 
caused by instability could increase abundances of these 
elements in the atmosphere and make the corresponding 
spectral features detectable. 

Turbulent motions can also lead to an efficient hori- 
zontal diffusion of the accreted material. The characteris- 
tic timescale of turbulent diffusion over the neutron star 
surface can be estimated as t s ~ L 2 /vtx- Substituting 
expression (67) and assuming e <~ H/L, we obtain 



t. 



£>13 



(69) 



If \ Z /X X = k x /k z <~ H/L and Bi 3 <~ 0.1, x <~ 1 then 
spreading of the accreted material over the surface pro- 
ceeds on a timescale of about a few hours. 

Apart from mixing, turbulence can also enhance the 
heat transport. The turbulent thermal diffusivity is com- 
parable to vtz-, and the ratio of turbulent and electron 
thermal diffusivities is 



X(0) 



B 13 T 7 - 1 (l + x 2 ^)x- 1 \^. 



(70) 
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In general, these quantities can have the same order of 
magnitude. Since VT is subadiabatic, turbulent motions 
increase the difference between the surface and internal 
temperature. Note that turbulent heat transport along the 
surface and vertically should be different for the consid- 
ered instability with more efficient diffusion along the sur- 
face, which can reduce the surface temperature gradient 
and decrease the contrast between the polar and equato- 
rial temperature. 

The non-oscillatory instability, despite being slower 
than the oscillatory one, can also manifest itself in accret- 
ing neutron stars. A sufficiently strong magnetic field can 
efficiently funnel the accreted material onto the magnetic 
pole. The characteristic radius of such accretion spots in 
the polar region is likely to be substantially smaller than 
the stellar radius, L <~ 10 4 — 10 5 cm. In a moderate mag- 
netic field satisfying the inequality Lx/Hx > a > 1, the 
growth rate can be estimated if we suppose in Eq. (62) 
that k x /k z ~ aHx/Lx (see Eq. (57)). Then, 

7f -lO^l + z 2 / 3 )- 1 ^-^ s- (71) 

If Z is not very large and the radius of accretion spots on 
the surface is ~ 100 m, then the growth time of instability 
can reach 10 4 — 10 5 s. In the non-linear regime, the tur- 
bulent velocity induced by this instability is rather small, 
~ 0.1 — O.OlL^g cm/s, if x ~ 514 ~ Hx3 ~ 1. Slow tur- 
bulent motions induced by the non-oscillatory instability 
in the spot can produce slow variations in the X-ray lu- 
minocity of accreting pulsars on the same timescale, which 
should be longer for accreting neutron stars with a weaker 
magnetic field. 

Note that the criteria and growth rates of instability 
have been derived under the assumption c s > ca or p > 
B 2 /4ir. Neglecting the quantum effects, we havep w 9.7 x 
10 22 a; 5 / 3 dyne/cm 2 in the region where the electron gas is 
non-relativistic and x < 1. Then, our results applies if 

B < B cr (x) ss 1.1 x 10 12 a; 5 / 6 G. (72) 

Stability properties of the layers where ca > c s (or 
B > B cr ) arc different from those considered in our study. 
Generally, perturbations of T and X can be suppressed 
for modes with 7 <~ lua > w s . As a result, a strong mag- 
netic field in the region where ca > c s can reduce the 
growth rate of instabilities by a factor (ca/c s ) 2 . Therefore, 
the turbulent transport on a short time-scale can be sup- 
pressed in these layers, and instability is less efficient in 
the surface region where ca > c s or 

p < p cr w 2.8 x 10 7 B^ 5 g/cm 3 . (73) 

A difference in turbulent transport between the regions 
with p > p cr and p < p cr could produce a qualitative 
difference in the burst activity of strongly and weakly 
magnetized accreting neutron stars. The thermonuclear 
X-ray burst is triggered when H/Hc burning becomes un- 
stable at the base of the accumulated layer, at a density 
<~ 10 5 - 10 6 g/cm 3 (see, e.g., Bildsten 1998, Cumming 



& Bildsten 2000). However, strongly magnetized accret- 
ing neutron stars are not observed among X-ray bursters, 
despite accretion at rates for which the burning of H/Hc 
should be unstable. This is usually considered as evidence 
that material is not spreading over the surface on short 
timescales, so the local accretion rate is high enough to 
stabilize the burning (Joss & Li 1980). In this model, the 
burning can be stabilized if spreading is suppressed only 
in the surface layer above the region with p ~ 10 5 — 10 6 
g/cm 3 (or x < 0.1 — 1) where the H/He burning occurs. 
One can estimate from Eq. (72) that instability and mix- 
ing become less efficient in the layers with p < 10 5 — 10 6 
g/cm 3 if 

B > B cr {x ~ 0.1 - 1) - 10 11 - 10 12 G. (74) 

Therefore, one can expect that the accreted material is 
spreading over the surface layer with p < 10 5 — 10 6 g/cm 3 
on a short timescale only in weakly magnetized neutron 
stars with B < 10 11 - 10 12 G, and the H/He burning in 
such stars should be unstable. On the contrary, in strongly 
magnetized neutron stars with B > 10 11 — 10 12 G, the 
spreading over the surface layer requires a longer time, 
and these neutron stars will not manifest themselves as X- 
ray bursters. The difference between strongly and weakly 
magnetized accreting neutron stars will be considered in 
more detail in a forthcoming paper. 
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